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In spatially distributed cellular systems, it is often convenient to represent complicated auxiliary 
pathways and spatial transport by time-delayed reaction rates. Furthermore, many of the reactants 
appear in low numbers necessitating a probabilistic description. The coupling of delayed rates with 
stochastic dynamics leads to a probability conservation equation characterizing a non-Markovian 
process. A systematic approximation is derived that incorporates the effect of delayed rates on 
the characterization of molecular noise, valid in the limit of long delay time. By way of a simple 
example, we show that delayed reaction dynamics can only increase intrinsic fluctuations about the 
steady-state. The method is general enough to accommodate nonlinear transition rates, allowing 
characterization of fluctuations around a delay-induced limit cycle. 



I. INTRODUCTION 

Biochemical circuits underlying many complicated cell 
functions, disease states, or viral propagation are often 
modeled by systems of delayed differential equations [TJ 
EJ O SJ [3 |6] , where the delay time represents auxiliary 
reaction pathways or spatial transport. It is well-known 
that in chemical reaction networks at the cellular level, 
many of the reactants are present in low copy number and 
therefore intrinsic noise is simply one of the operating 
constraints 7 . It is, however, difficult to develop and 
analyze models that contain both stochastic dynamics 
and delayed reaction rates. 

For systems without delay, a stochastic description 
typically takes the form of a chemical Master equation 
that governs the probability distribution P(n, t) of find- 
ing the system in state n at time t conditioned upon 
some initial state [H [S] ■ The Master equation can rarely 
be solved exactly, and various methods have been devel- 
oped to approximate the evolution of P(n, t). Perhaps 
the most well-known is the Kramers-Moyal expansion, 
which when truncated after the second term results in a 
diffusion equation with nonlinear drift and non-constant 
diffusion, called the Fokker-Planck equation [TU], [TTJ . If 
the deterministic system evolves along a stable trajectory 
near a stable fixed point, van Kampen [12J has developed 
an alternate approximation of the Master equation that 
relics upon a perturbation expansion in some extensive 
quantity, providing a consistent characterization of the 
fluctuations in terms of a Fokker-Planck equation with 
linear drift and constant diffusion from whence the mean 
and variance are easily computed. (For a more detailed 
discussion, see [HUES] and references therein.) 

Nevertheless, in many systems the individual reaction 
events depend upon the past state of the system [71 H5]. 
and the methods developed to approximate the solution 
of the Master equation are no longer appropriate. Writ- 
ing the transition probability of moving from state n' 
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to state n in an interval dt as W n \ n dt and the two- 
point joint probability distribution of finding the system 
in state n at time t and in state m at time t — r as 
p2(n, t] m i t — t), the delayed dynamics introduce a con- 
volution term into the probability conservation equation, 

^Qt ^ = E Wn'.niV, t) ~ W n , n ,P(n, t)+ (1) 

11 7 

E E ^',n P 2 ( n '» ^ m > * - T ) ~ W n,„' ^2 (n, t; m, t - T ) . 
n' m 

Here W^, n (m) depends upon the state at a time r in the 
past: m = n(t — r). Eq. [I] is no longer a closed equation 
for P(n, t) since it now includes the unknown distribu- 
tion P2(n, t; m, t — r). As a consequence, it no longer 
describes a Markov process and methods used to treat 
the standard chemical Master equation require modifica- 
tion. 

Many investigations using stochastic simulation algo- 
rithms have illustrated the importance of intrinsic noise 
in systems with delay [3J |SJ [7] , while past analytic work 
has focused primarily upon approximations of the de- 
layed nonlinear Fokker-Planck equation [UMTS], stochas- 
tic delay differential equations [5J [TSJ [T7] or exactly- 
solvable random-walk models [TSJ \T§\ . While each pro- 
vides considerable insight into the interdependent effects 
of noise and time delay, comparatively little work has 
been done to connect the underlying discrete probabil- 
ity conservation equation to these continuous approxi- 
mations. In what follows a perturbation scheme is de- 
veloped that, under the condition that the delay time 
exceeds the relaxation time of the deterministic system, 
allows a general probability conservation equation to be 
approximated by a delayed linear Fokker-Planck equa- 
tion, thereby making connection to past studies. Recent 
work by Bratsun et ai. |13j has explored very similar 
questions, although their method is applicable only to 
systems with linear reaction rates. The delayed linear 
noise approximation, which is an extension of van Kam- 
pen's linear noise approximation [12 , provides a consis- 
tent characterization of intrinsic fluctuations in delayed 
systems and is sufficient to show that under fairly general 
conditions delayed reaction events can only increase the 
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magnitude of these fluctuations. A simple example of a 
birth/death process is used to provide a concrete imple- 
mentation of the method, and a nonlinear predator-prey 
model illustrates characterization of fluctuations along a 
delay-induced limit cycle. 

II. MATHEMATICAL METHODS 

A stochastic model of a network of chemical reactions 
governs the probability distribution P(n, t) of finding the 
system in state n at time t, with dynamics given in terms 
of the stoichiometric change resulting from the comple- 
tion of each reaction and the propensity of occurrence for 
each reaction event, recorded, respectively, in the stoi- 
chiometry matrix S and the propensity vector u |201 121j . 



Consider a system with N reactants that can combine 
through one of R reactions. To facilitate the inclusion of 
delayed kinetics into the formalism, we separate the R 
reactions into two groups: those with rates {i>j (n)}j^ 1 
depending upon the current state of the system and those 
with rates {evj (m, r\)}^ =R +1 depending upon the past 
state of the system m = n [t — r), where r is the delay 
time. The parameter e measures the delayed feedback 
strength; throughout, we shall assume that the feedback 
is weak, and explicitly retain leading-order terms in e. 

To keep the notation compact, we introduce the step- 

—S' ■ 

operator E ; ,J that acts to increment the variable rii by 

a. . 

an integer — SV,-: E i — f( n i~ Sij)- Denoting the 

system volume by V, Eq.[T] takes the form [T31 |2"T1 12"2"] . 
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delayed - state dynamics 
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where ({rii}) is the Heaviside-step function that en- 
sures no delayed reaction occurs if completion results in 
the unphysical end-state < for any elements of n. 
Throughout, only non-consuming reactions are consid- 
ered O [23] , i.e. reactants of an unfinished reaction are 
allowed to participate in new reactions. 

The solution of the full distribution P (n, t) is not pos- 
sible in general, therefore we seek an approximate solu- 
tion. To that end, we make the ansatz that the number 
of react ant molecules is large enough that the discrete 
molecule numbers ??, can be represented by the continu- 
ous deterministic concentrations Xi and some continuous 
fluctuations on that scale as the square-root of the num- 
ber of molecules [12] , 

tl{ = Vxi + W a>i and mj = Vx\ + W fy, (3) 

where V is the system volume and x\ = Xi(t — t). Us- 
ing the auxiliary variable fli, the delayed fluctuations are 
written as fli (t — r) = (t — t) to emphasize the ap- 
proximation made below; specifically, that the delay time 
is sufficiently large that at and (3i can be treated as inde- 
pendent random functions. As we show below (Eq. 18 1, 
that assumption is consistent with the requirement that 
the delay time is much larger than the characteristic re- 
laxation time of the deterministic equations and that the 
delayed feedback is weak (e « 1). In the expansion be- 
low, 1/V is assumed small, although since Xi is held fixed, 
(equivalcntly, one assumes that ni is large). The result- 
ing approximation will be called the delayed linear noise 
approximation. 

Invoking the linear noise approximation by Taylor- 
expanding the microscopic transition rates about the 



macroscopic trajectory x (t) in powers of \/V [12 |2"T] . 
we have 



iV 



Uj (m, n, V) w Vj (x T , x) + -= ^ 



dxi 



dv 

dx 



with an analogous expression for Vj (n, V). The rates Vj 
correspond to the deterministic reaction rates (see [21] 
for a discussion of the difference between and v). The 
discrete step-operator E is likewise expressed as a Taylor 

series in VV' 1 [H[5T], 

N I N 1 N 



i=l 



,k=l 



where di — d/dcti. The one-point P(n,t) and two- 
point joint probability P2 (n, t; m, t — r) can be written 
in terms of the single distribution and joint distribu- 
tion of the fluctuations about the macroscopic trajectory, 
II (a, t) and n 2 (a, t; f3,t — r), respectively, via the linear 
change of variables suggested by Eq. [3] 

P(n,t) = V- N ^ 2 U( a ,t), (4) 
P 2 (n, t; m, t - r) = V- N U 2 (a, t; /3, t - r) , (5) 

where a and (3 are centered upon x (t) and x (t — r) , 
respectively, and the factor involving V comes from the 
normalization of II, 



/oo />oo 
P(n,t)dn=V- N/2 Tl(a ) t)da=l. 
-OO J —OO 



(6) 



3 



The decoupling of the deterministic component x from 
the stochastic component a, with a \j\/V scaling of the 
fluctuations implied by the ansatz Eq. [3] is the most 
fundamental step in the approximation. That assump- 
tion leads directly to the normalization above, and al- 
lows subsequent terms in the perturbation expansion to 
be ordered in terms of powers of X/yV [IS]. 

For long delay time t, and small delayed- feedback 
strength s, the fluctuations at time t are approximately 
independent of the fluctuations at time t — r, allowing 
the joint-distribution to be factored as, 

n 2 (a, t; (3, T) rs II (a, t) x II (/3, T) + eh(a, P) (7) 
T < t - T, 

where h(at,f3) must obey the consistency condition, 
J h(a,(3)da = J h(a, f3)d(3 = 0. 



Notice the factoring of the fluctuations II2 is not equiv- 
alent with assuming independence in the full state, 
P(n, t;m, t) = P(n,t) x P(m, r) (as is done in [13]). 
Independence of the full state is inconsistent with the 
deterministic rate equations, 



dx 

~dt 



S • v = f (x r ,x) , 



(8) 



where, in the limit of large n, the present state is com- 
pletely determined by the past-states (except, perhaps, at 
steady-state). Moreover, independence of the full state 
implies that higher-order correlations, including the au- 
tocorrelation function K(t — s), vanish for \t — s\ > t. 
We show below (Eq. 18 ) that this is not the case. 

In the limit 1/V — ► (n, — > 00), the Heaviside step 
function is = 1; and with the factored joint 

distribution, Eq. |7J the sum over m is replaced by the 
integral over f3, 
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N r 



dx 
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dv 
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Ha: 



y-n(a,f) 



Here, the first two terms on the right-hand side follow 
from the normalization condition on TI(/3, t — r), Eq. |6] 
and the third from, 

/■oo 

v- N ' 2 / ftn(p,t-T)dp = {Pi). 



Substituting the expanded terms in Eq. [2] using the 
chain-rule to write the partial derivatives of P(n, t) in 
terms of II and a. [12] , and collecting in powers of VV, 
the zero'th order term is simply the deterministic de- 
layed reaction rate equations, Eq.[sJ At \[V , we obtain 
the equation characterizing the probability distribution 
of the fluctuations a(t), 

dU 

~dt 



-T iS di (n ; ll) 



Di 



d^U 



£ r-.(^)^n, (9) 



where, 



r, 



d[S ■ v]i 



, C-L l? 



d[S-v]i 



d Xj ' 13 dx] 



D = S diag [u] ■ S 1 



Eq.[9]is a closed diffusion equation for II(a, t) with coeffi- 
cients that are linear in the fluctuation variables a. The 
matrix D represents the diffusive effects of the fluctua- 
tions, while the matrices T and T T represent the restora- 
tive drift in the system [22] |24] . Eq. [9] is not quite a 
Fokker-Planck equation since it contains the delayed av- 
erage (j3j) = (ctj(t — t)) in the drift coefficient. Never- 
theless, the initial conditions can be chosen so that the 



last term in Eq. [9] vanishes, as we now show. Multiplying 
^ by en and integrating yields the evolution equation 
for the mean. 



<*(«(*)) 
dt 



T • (a (t))+eT T ■ (a(t-r)) 



(10) 



The initial average (ot(t)) can always be absorbed into 
the initial conditions on x (t) so that (a (t)) = for t ^ 0, 
thereby ensuring that (a (t)) — for all time. Without 
loss of generality, then, we write Eq. [9] as the Fokker- 
Planck equation with coefficients linear in a, 



dt 



= -J2 FiA K-n) + i^A'hu. (11) 



It is important to note that although the fluctuations at 
time t are independent of fluctuations in the past, they 
are conditioned by the macroscopic solution x (t) (and 
x (t — t)) through the coefficient matrices T and D. 

A consequence of Eq. llXl is that, to O (V -1 ) , the 
fluctuations are Gaussian distributed with covariancc 
OiiOij) — (oti) (oij) = (aiOij). Multiplying Eq. 11 by 



cuOLj and integrating over all a. yields a dynamic equation 
governing S [22, p. 211], 



dt 



r-H + s-rt + D, 



(12) 



(where is the matrix transpose of T; not to be con- 
fused with T T ). At steady-state, the coefficient matrices 
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r and D (and therefore S) will be constant, satisfying 
the fluctuation-dissipation relation, 



Laplace transform of the mean, 



r-H + s-r t + D = o. 



(13) 



The diffusion matrix D is symmetric and positive semi- 
definite by construction, so that a steady-state probabil- 
ity distribution is only possible if the drift term T bal- 
ances the diffusion D. With long-delay in the reaction 
kinetics, the restorative influence of T T no longer appears 
in the equation governing the fluctuations (Eq. Ill, so 



that although the delayed dynamics increase the magni- 
tude of the diffusion matrix _D, the dissipation due to 
r T is lost. Therefore, in the limit of long delay time, de- 
layed dynamics can only serve to increase the magnitude 
of intrinsic fluctuations. 



III. STEADY-STATE AUTOCORRELATION 
FUNCTION AND SPECTRUM 

The fluctuation-dissipation relation, Eq. |13[ and the 
evolution equation for the mean, Eq. |10[ together pro- 
vide an expression for the time-autocorrelation matrix 
for the fluctuations about the steady-state, K(t) = (a(t)- 
ct T (0)). The steady-state autocorrelation function K(i) 
is, by definition, 



«(*)>a(0) -a T (0)n(a,0)rfa, (14) 



K (i) = jja.' (t) ■ a T (0) n 2 (a', t; a, 0) da! da. 
Re-writing in terms of the conditional probability, 
K (t) = J J a' (t) ■ a T (0) n (a', t\a, 0) II (a, 0) da' da 
K(t) = 

where (a (t)) a , \ is the solution of Eq. 10 with initial 
condition a(0), and II is the equilibrium distribution of 
a (o) mug. Eq.m requires that the conditional prob- 
ability density H(a\t\a,0) be written as a perturbation 
expansion in e, 

II(a', t\a, 0) = n°(a', t\a, 0) + ell^a:', t\a, 0) + 0(e 2 ). 

(15) 

Consequently, in the conditional average (a (t)) a ^, only 
terms to first-order in e are retained. 

The conditional average (a (t)) a ^ is obtained from 



Eq. 10 which is easily solved via Laplace transform. The 
equilibrium correlation function is an even function of 
time-difference alone, equivalent to boundary condition 
(a (t)) =0 for t < 0, leading to the formal solution 



(&00) t 



(0) 



[si - r 



ee 



(a(0)). (16) 



The derivation of Eq. [TTJ assumes e — > 0, so to remain 
consistent, we retain only leading-order terms in e in the 



(a( S )) = [ S I-T}- 1 -(a(0)) 



-ee 



si-ri 



r T • [si - n 



(a(0)>. (17) 



With substitution into Eq. [TT] using the fluctuation- 
dissipation relation, Eq. |13| the autocorrelation function 
is 

K(t) = exp [Ft] ■ 3 S (18) 

+ee(t - t)£-} t {[si - rr 1 ■ v T ■ [si - rr 1 } ■ s„ 

where Q(t — r) is the Heaviside step function, and £^~_ T 
is the time-shifted inverse Laplace transform. The first 
term produces an exponential drop from t — 0, while 
the second term produces an anti-correlated second peak 
slightly beyond t = r; higher-order terms in e produce al- 
ternating correlated/ anti-correlated peaks of magnitude 
0(e n ) for the n th peak. 

The fluctuation spectrum follows immediately from the 
autocorrelation function. We denote the £-correction to 
the autocorrelation function K corr (s), 

k corr ( s ) = e - ST [si - rr 1 • r r • [si - rr 1 , (19) 

then the spectrum S(cj) is [TT] . 

poo p0 

S(u)= e~ iut K(t)dt+ e iwt K^(-t)dt 

JO J -00 

= H^i + rp 1 d [iuji + r^y 1 

+e {K corr (icj) ■ E s + E s ■ Kl orr (-iw)} , (20) 

where D is the diffusion matrix introduced in Eq. [^eval- 
uated at the steady-state. 



IV. LINEAR EXAMPLE - DELAYED PROTEIN 
DEGRADATION 

To provide some context for the formal derivation 
above, consider a simple birth/death model with delayed 
degradation |13| . The total number of species X evolves 
via the following three reactions, 



constant synthesis: X — > X + 1 
linear degradation: X -4 X — 1 
delayed degradation: X =4 X — 1 



v\ = 1, 
v-x = 5 ■ x, 
z^3 — eS ■ x T . 



(21) 



The reaction rate vector (in units of concentra- 
tion/time) is given by v = [7, S • x, eS ■ x T ] and the sto- 
ichiometry matrix is S = [1,-1,-1]. The determin- 
istic reaction rate equation for the concentration x(t) 
is then governed by the delayed differential equation, 
x = S • v = 7 — 5 • x — sS ■ x T . The auxiliary coef- 
ficient matrices in Eq. [TT] are the scalars T = —6 and 
D = 7 + 5 ■ x + eS ■ x T . For the sake of simplicity, we con- 
sider the fluctuations about the steady-state x s , where 
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FIG. 1: Fano factor as a function of delay time. The 

t — (dotted) and r — > oo (solid) approximations to the Fano 



factor, var[a;]/(a;} (Eq. 23 1, provide lower and upper bounds 
on the estimate of intrinsic noise over the range of delay times, 
as compared to stochastic simulation data generated from the 
model shown in Eq. 21 (filled circles). The delayed linear noise 



approximation adequately characterizes the fluctuations even 
when the delay time is of the order of the other time scales 
in the problem (i.e. tS « 1). Here, the deterministic reaction 
rates are 7 = 100, S = 4, e — 1/4, with V = 1. Simulation 
data is from an ensemble of 10° members using the delayed 
Gillespie algorithm [23] . 



x s = 7/(5(1 + e). From Eq. 13 the variance of the fluc- 
tuations about x s is, (using Eq. [3]), 



- (x.r = 



V 



V 



V 8' 



A useful measure of the relative magnitude of the fluctu- 
ations is the fractional deviation rj 2 ., 



V -xl 



1 



(1 + e), 



(22) 



where here N s — V ■ x s is the number of molecules in 
the steady-state. Without delay, this simple example de- 
scribes a Poisson process, and for a Poisson process the 
Fano factor F — rj 2 N s = E/x s = 1. From Eq. [22 



F = 



1, T 

1 + e, t- 



•> 
00. 



(23) 



The Fano factor is a particularly convenient statistic to 
contrast the ordinary and delayed linear noise approx- 
imations, as well as illustrating the delay time neces- 
sary for the present approximation to hold. Figure [l] 
shows the Fano factor F estimated from stochastic sim- 
ulation [201 HI] with 7 = 100, 8 = 4 and e = 0.25, as 
compared with the long delay time (solid) and short de- 
lay time (dotted) estimates. Notice the cross-over occurs 
for r w that is for delay time comparable to the 
natural time scale of the undelayed kinetics. 



The autocorrelation function is given by Eq. [18] 



K(t) = e~ 5t {l - eB (t- r) Se Sr (t - r)} 



7 
6' 



(24) 



This expression coincides with the result of Bratsun et 
al. [T3] to 0(e), and as they demonstrate, K(t) very 
faithfully reproduces the autocorrelation from simulation 
data. Furthermore, the autocorrelation can be used to 
identify 'quasi-cyles' where regular oscillations emerge 
from deterministically stable systems [T7J ISZl HE] ■ 

The delayed-degradation model is used as a transpar- 
ent illustration of the method, but the same results can 
be obtained by other methods (for example, via moment- 
generating functions [13]). In contrast, the model in the 
next section contains more realistic, nonlinear transition 
rates, and consequently cannot be treated by existing 
methods. Yet nonlinear rates abound in physical applica- 
tion and exhibit rich dynamics, as the following example 
demonstrates. 



V. NONLINEAR EXAMPLE - 
PREDATOR-PREY DYNAMICS 

The methodology outlined in Section [n] makes no as- 
sumptions about the nonlinearity of the transition prob- 
abilities in the stochastic model, opening up the possibil- 
ity to study the dynamics of delayed nonlinear systems. 
As an example capable of exhibiting asymptotic stabil- 
ity and limit cycle behavior, consider the predator-prey 
model with delayed predator birth, 



dP 
~dt 



a ■ P(t) 



K 



P(t) 2 -b-P(t)Z(t), 



dZ 
~dt 



C- P(t-T)Z(t-T) -d-Z(t), 



where P(t) is the density of prey, Z(t) is the density of 
predators and K is the carrying capacity of the environ- 
ment. Here, r is the delay time associated with gestation 
before the birth of predators. Assuming each birth event 
produces a litter of 1, the reaction network, in volume V, 
takes the form. 



(25) 
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P + 


1; 


v x = a ■ 2£ 
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P- 


1; 


a n p np — 1 
2 K V V 


p 


P- 


1; 




Z 


Vj 


z- 


1; 


u 4 = d-^ 


z 




z + 


1; 


I*-C- np( y T) ■ nZ( y T) 



where n P /V = P(t) and n z /V = Z(t). 
With a suitable nondimensionalization, 

to = l/a, Pq = K , Zq = a/b, e = c ■ K/a, and 8 = d/a, 

the deterministic model equations reduce to, 

dP 

— = P(t)~P(t) 2 -P(t)Z(t), 
dZ 

— = e-P(t~T)Z(t-T)-8-Z(t), (26) 
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A. Asymptotically stable 

For (e = 0.15,5 = 0.05), the system remains asymp- 
totically stable in both limits, r = (Figure [2]A_) and 
r = 30 (Figure [2(3). From Eq. 13 long delay time re- 



7.0 7.5 
P(x10 3 \ 



duces the stability imparted to the system through T. 
As a consequence, the variance of the fluctuations is ex- 
pected to increase with increasing delay time r. This 
increase is evident in the stationary probability distribu- 
tion of the fluctuations derived from stochastic simula- 
tion (Figure pj. The ellipses shown in the figure corre- 
spond to the first- and second-standard deviations of the 
steady-state Gaussian distribution predicted by the ap- 
proximation, while the density plot represents extensive 
stochastic simulation data generated using Gillespie's al- 
gorithm [mug. 

The most striking consequence of the delay on the 
intrinsic fluctuations is the increased magnitude of the 
cross-correlation between P and Z. As r — * 00, the de- 
layed rate c-P(t—T)Z(t—T) no longer offers compensation 
to the predation event with rate — b ■ P(t)Z(t), resulting 
in the increased cross-correlation. This is an example 
of the nontrivial effect of delayed dynamics on intrinsic 
fluctuations, even though the equilibrium point is sta- 
ble. In situations where the delay affects not only the 
fluctuations, but the underlying stability itself (as is the 
case for delay induced limit cycles) , the analysis becomes 
more complicated. 



B. Limit cycle 



FIG. 2: Steady-state fluctuations in a nonlinear model 

(e = 0.15, S = 0.05). A. Density plot of the equilibrium fluc- 
tuations from stochastic simulation (10 6 realizations). Solid 
curves correspond to the first- and second-standard deviation 
ellipse computed by an ordinary application of the linear noise 
approximation [TJ]. B. Same model parameters as in panel 
A, but with delayed predator birth (r = 30). The solid curves 
correspond to the first- and second-standard deviation ellipse 
computed by the delayed linear noise approximation. 



where r has been likewise nondimensionalized by l/a. 
The equilibrium point corresponding to coexistence of 



the populations is (P, Z) 



' , 1 — -), leading to a nec- 



essary condition for stable coexistence, with and without 
delay, that S < e. 

It is well-known that delayed rates can have a desta- 
bilizing effect on population dynamics |29j . and in fact 
can generate limit-cycles in otherwise stable models [6, 
I30"| I3T] . To illustrate the approximation method and the 
destabilizing effects of delay, we consider two values for 
the delay time, t = and r = 30, in two parameter 
regimes - the first chosen so that the equilibrium remains 
asymptotically stable for both values of the delay time, 
the second chosen so that a limit-cycle appears for large 
delay r. 



For (e = 0.25,(5 = 0.05), the system is asymptotically 
stable for r = 0, but a limit-cycle appears for r = 30. 
By separating the fluctuations tangent to the limit cycle 
from those transverse [UJ [331 [32 [33], the delayed lin- 
ear noise approximation is easily extended to a system 
exhibiting a limit-cycle. 

Briefly, a moving coordinate frame is introduced using 
as a basis the unit vectors tangent (s) and normal (f ) to 
the limit cycle. In the moving frame, the covariance of 
the transverse fluctuations S rr decouples from the diver- 
gent fluctuations along the limit cycle, and is character- 
ized by a stable evolution equation (cf. Eq. 12), 



dt 



or' s _i_ n' 

rr ' — Tr i ^rr' 



(27) 



where T' rr and D' rr are elements of the drift and diffusion 
matrices in the moving frame, 

r^U-r-Tjt-r- U f , D' = U D U f , (28) 

dt 

and U is the rotation matrix generated from the deter- 
ministic rate equations f(x,x T ), 



U = 



1 



/l — J 2 
h h 



(29) 
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FIG. 3: Fluctuations around a delay-induced limit cy- 
cle (e = 0.25,5 = 0.05). A. Envelope of the standard- 
deviation of the fluctuations transverse to the limit cycle, as 
computed by Eq|27| using the delayed linear noise approxima- 
tion. B. Stochastic simulation of the system along the limit 
cycle. 



Figure [3] illustrates the estimate of the fluctuations 
along the delay-induced limit cycle via the delayed lin- 
ear noise approximation (Figure ) , compared with the 
result of a stochastic simulation (Figure [3)3 )■ The width 
of the envelope of the fluctuations is not uniform around 
the orbit, reflecting the state-dependent drift T and dif- 
fusion D matrices in Eq. |12| This same non-uniformity 
is also observed in the stochastic trajectory. 

The nonlinear predator-prey model demonstrates the 
utility and comparative simplicity of the delayed linear 
noise approximation - once the network is written in 
terms of the stoichiometry matrix and the propensity 



vector, despite the lengthy derivation, Eqs. [9] [12] and 
18 allow algorithmic characterization of the fluctuations. 
VI. DISCUSSION 



In models of cellular chemical reaction systems, spatial 
transport and long auxiliary pathways are often repre- 
sented using time-delayed reaction rates. At the meso- 
scopic level, delayed dynamics result in a probability con- 
servation equation that characterizes a non-Markovian 
process. Since analytic solutions are rare, approximation 
of the governing equations are necessary. In the limit of 
large numbers of molecules, weak delayed feedback and 
long delay time, we have derived the leading order behav- 
ior of a probability conservation equation with delayed 
transition rates from an expansion in the system volume 
V. The fluctuations are characterized by a linear Fokker- 
Planck equation, in accordance with the linear noise ap- 
proximation of the undelayed case |12j . and coinciding 
with a delayed random walk in a quadratic potential |19j . 
We find that the delayed dynamics contribute unevenly 
to the drift and diffusion coefficients of the Fokker-Planck 
equation, and conclude that long time-delay can only in- 
crease the magnitude of intrinsic fluctuations for systems 
where the delayed feedback has a stabilizing effect. 

Here, we have focused upon two example systems - one 
that evolves toward a stable steady-state, the second is 
a nonlinear model exhibiting a delay-induced limit cycle. 
It is often the case that models with delayed rates are 
used to describe oscillatory dynamics [3]. The delayed 
linear noise approximation is easily adapted to systems 
evolving along a stable limit cycle by a simple change of 
coordinates [2"H 155] . 

Finally, the effect of noise on the macroscopic behavior 
of a system is not always additive, and in fact noise can 
generate ordered oscillations from a deterministically sta- 
ble model [HI [34] . These noise induced oscillations have 
been proposed as a mechanism to extend the capacity of a 
given network to sustain oscillations [35 , 36 . The results 
derived above, specifically the autocorrelation function 
Eq. [18] allow the method developed for studying noise- 
induced oscillations in undelayed systems to be applied 
to systems characterized by delayed dynamics [25] . 
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